/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

shornuf : Shor's factoring alorithm implemented using the qdd library

Version 0.1 

Copyright (C) 1999 David Greve <dagreve@iname.com>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

#include "qdd.h"

// DAG -- something wrong here ..
extern void qdd_init(int,int);
extern Integer Ipow(long, long);


/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  shornuf : Shor's factoring algorithm implemented using the qdd
  library.  This program is translated almost line for line from
  Bernhard Oemer's "QCL" implementation of this same algorithm.

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

void lt (QuantumRegister & B, const Integer & A, QuantumBit & FLAG, QuantumRegister & CONT, int dir) {

  int size = B.Size;

  if (dir > 0) {

    ~ FLAG;                        // Start by assuming (B<A) .. set flag = true.

    int i = size - 1;

    if (BIT(A,i)) {
      CONT[i-1] ^=    B[i]  ;      // If A[i] is true, continue only if B[i] is also true (tie).
    }
    else {
      CONT[i-1] ^= (! B[i]) ;      // If A[i] is false, continue only if B[i] is also false.
      FLAG      ^=    B[i]  ;      // also, If B[i] is true, B cannot be less than A.
    }

    for (i=(size - 2);i>=1;i--) {
      if (BIT(A,i)) {
	CONT[i-1] ^= CONT[i] &    B[i]  ;    // If A[i] is true, continue only if B[i] is also true (tie).
      }
      else {
	CONT[i-1] ^= CONT[i] & (! B[i]) ;    // If A[i] is false, continue only if B[i] is also false.
	FLAG      ^= CONT[i] &    B[i]  ;    // also, If B[i] is true, B cannot be less than A.
      }
    }
    
    if (BIT(A,i)) {                // OK, last bit ..
      FLAG ^= CONT[i] &   B[i]  ; // If A[i] is true and we are still checking, B < A iff B[i] = 0.
    }
    else {
      FLAG ^= CONT[i] ;            // If A[i] is false, B cannot be less than A if we are still checking.
    }

  }

  else { // reverse

    int i = 0;

    if (BIT(A,i)) {                
      FLAG ^= CONT[i] &   B[i]  ; 
    }
    else {
      FLAG ^= CONT[i] ;           
    }

    for (i=1;i<=(size - 2);i++) {
      if (BIT(A,i)) {
	CONT[i-1] ^= CONT[i] &    B[i]  ;    
      }		               
      else {	               
	CONT[i-1] ^= CONT[i] & (! B[i]) ;    
	FLAG      ^= CONT[i] &    B[i]  ;    
      }
    }
    
    if (BIT(A,i)) {
      CONT[i-1] ^=    B[i]  ;
    }
    else {
      CONT[i-1] ^= (! B[i]) ;
      FLAG      ^=    B[i]  ;
    }

    ~ FLAG;

  }

}

void qcl_muxhaddbit(   int a0,    int a1, QuantumBit & sel, QuantumBit & b, QuantumBit & ci_sum, QuantumBit & e, int dir) {

  if (dir > 0) {

    if (a0 != a1) {           // a0 and a1 differ?
      ci_sum ^= (! sel) & e;  // add a
    } 
    else {
      if (a0 && a1) { 
	ci_sum ^= e;         // add a
      }
    }

    ci_sum ^= b;             // add b

  }

  else {

    ci_sum ^= b;             // add b

    if (a0 != a1) {           // a0 and a1 differ?
      ci_sum ^= (! sel) & e;  // add a
    } 
    else {
      if (a0 && a1) { 
	ci_sum ^= e;         // add a
      }
    }

  }

}

void muxhaddbit (int a0, int a1, QuantumBit & sel, QuantumBit & b, QuantumBit & ci_sum, QuantumBit & e, int dir) {

  if (dir > 0) {

    QuantumBit M(0);

    if (a1) M ^= e & sel;

    if (a0) M ^= e & (! sel);

    ci_sum ^= e & M;
    ci_sum ^= e & b;

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

    try {
      M.unMixed();
    }
    catch (QuantumException e) {
      cerr << "muxhaddbit() : Mixed Bits\n" ;
      throw e;
    }

  }

  else {

    QuantumBit M(0);

    if (a1) M ^= e & sel;

    if (a0) M ^= e & (! sel);

    ci_sum ^= e & b;
    ci_sum ^= e & M;

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

    try {
      M.unMixed();
    }
    catch (QuantumException e) {
      cerr << "muxhaddbit() : Mixed Bits\n" ;
      throw e;
    }

  }

}

void muxaddbit (int a0, int a1, QuantumBit & sel, QuantumBit & b, QuantumBit & ci_sum, QuantumBit & co, QuantumBit & e, int dir) {

  if (dir > 0) {

    QuantumBit M(0);

    if (a1) M ^= e & sel;

    if (a0) M ^= e & (! sel);

    ci_sum ^= e & M;
    ci_sum ^= e & b;

    co  ^= e & M & (! ci_sum);
    co  ^= e & (! M) & b & (! ci_sum);
    co  ^= e & M & b & ci_sum ;

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

    try {
      M.unMixed();
    }
    catch (QuantumException e) {
      cerr << "muxaddbit() : Mixed Bits\n" ;
      throw e;
    }

  }

  else {

    QuantumBit M(0);

    if (a1) M ^= e & sel;

    if (a0) M ^= e & (! sel);

    co  ^= e & M & b & ci_sum ;
    co  ^= e & (! M) & b & (! ci_sum);
    co  ^= e & M & (! ci_sum);

    ci_sum ^= e & b;
    ci_sum ^= e & M;

    if (a0) M ^= e & (! sel);

    if (a1) M ^= e & sel;

    try {
      M.unMixed();
    }
    catch (QuantumException e) {
      cerr << "muxaddbit(-1) : Mixed Bits\n" ;
      throw e;
    }

  }

}

void muxaddreg (const Integer & a0, const Integer & a1, QuantumBit & sel, QuantumRegister & b, QuantumRegister & sum, QuantumBit & e, int dir) {

  qdd_assert(b.Size <= sum.Size, "Sum register too small");

  int size = b.Size;

  if (dir > 0) {

    int i;
    for (i=0;i<=size-2;i++) {
      muxaddbit(BIT(a0,i),BIT(a1,i),sel,b[i],sum[i],sum[i+1],e, 1);
    }

    muxhaddbit(BIT(a0,i),BIT(a1,i),sel,b[i],sum[i],e,1);

  }

  else {

    int i = size-1;

    muxhaddbit(BIT(a0,i),BIT(a1,i),sel,b[i],sum[i],e,1);

    for (i=size-2;i>=0;i--) {
      muxaddbit(BIT(a0,i),BIT(a1,i),sel,b[i],sum[i],sum[i+1],e,-1);
    }

  }

}

void addn (const Integer & a, const Integer & n, QuantumRegister & b, QuantumBit & flag, QuantumRegister & sum, QuantumBit & e, int dir) {

  int SumMaxBit = sum.Size - 1;

  if (dir > 0) {

    lt(b,n-a,sum[SumMaxBit],sum, 1);

    flag ^= e & sum[SumMaxBit] ;

    lt(b,n-a,sum[SumMaxBit],sum,-1);

    muxaddreg(Ipow(2,b.Size)+a-n,a,flag,b,sum,e, 1);

  }

  else {

    muxaddreg(Ipow(2,b.Size)+a-n,a,flag,b,sum,e,-1);

    lt(b,n-a,sum[SumMaxBit],sum,1);

    flag ^= e & sum[SumMaxBit] ;

    lt(b,n-a,sum[SumMaxBit],sum,-1);

  }

}

void swap (QuantumRegister & a, QuantumRegister & b) {

  a ^= b;
  b ^= a;
  a ^= b;

}

void eswap (QuantumRegister & a, QuantumRegister & b, QuantumBit & e) {

  int i;
  for (i=0;i<a.Size;i++) {
    a[i] ^= b[i] & e;
    b[i] ^= a[i] & e;
    a[i] ^= b[i] & e;
  }

}

void cxor (QuantumRegister & a, QuantumRegister & b, QuantumBit & e) {

  int i;
  for (i=0;i<a.Size;i++) {
    a[i] ^= b[i] & e;
  }

}

void oaddn (const Integer & a, const Integer & n, QuantumRegister & sum, QuantumBit & e, int dir) {

  QuantumRegister j(sum.Size);
  QuantumBit      f(0);

  if (dir > 0) {

    addn(  a,n,sum,f,j,e, 1);
    eswap(sum,j,e);
    f ^= e;
    addn(n-a,n,sum,f,j,e,-1);

    try {
      j.unMixed();
      f.unMixed();
    }
    catch (QuantumException e) {
      cerr << "oaddn() : Mixed Bits\n" ;
      throw e;
    }

  }
  else {

    addn(n-a,n,sum,f,j,e, 1);
    f ^= e;
    eswap(j,sum,e);
    addn(  a,n,sum,f,j,e,-1);

    try {
      j.unMixed();
      f.unMixed();
    }
    catch (QuantumException e) {
      cerr << "oaddn(-1) : Mixed Bits\n" ;
      throw e;
    }

  }

}

void muln(const Integer & a, const Integer & n, QuantumRegister & b, QuantumRegister & prod, QuantumBit & e, int dir) {

  QuantumBit T(0);

  if (dir > 0) {

    int i;
    for (i=0;i<=(prod.Size-1);i++) {
      if (BIT(a,i)) {
	prod[i] ^= b[0] & e;
      }
    }
    for (i=1;i<=(b.Size-1);i++) {
      T ^= b[i] & e;
      oaddn(((a << i) % n),n,prod,T, 1);
      T ^= b[i] & e;
    }

  }

  else {

    int i;
    for (i=(b.Size-1);i>=1;i--) {
      T ^= b[i] & e;
      oaddn(((a << i) % n),n,prod,T,-1);
      T ^= b[i] & e;
    }
    for (i=(prod.Size-1);i>=0;i--) {
      if (BIT(a,i)) {
	prod[i] ^= b[0] & e;
      }
    }

  }

}

void omuln(const Integer & a, const Integer & n, QuantumRegister & b, QuantumBit & e, int dir) {

  QuantumRegister j(b.Size);

  qdd_assert(((gcd(a,n) == 1) || (sign(a) == 0)), "Not Invertable");

  if (dir > 0) {

    muln(          a, n, b, j, e,  1);
    muln(invmod(a,n), n, j, b, e, -1);
    cxor(b,j,e);
    cxor(j,b,e);

  }
  else {

    cxor(j,b,e);
    cxor(b,j,e);
    muln(invmod(a,n), n, j, b, e,  1);
    muln(          a, n, b, j, e, -1);

  }

}

void expn(const Integer & a, const Integer & n, QuantumRegister & b, QuantumRegister & ex, int dir) {

  if (dir > 0) {

    ~ ex[0];

    int i;
    for (i=0;i<=b.Size-1;i++) {
      cout << "Bit : " << i << "\n" ;
      omuln(bitpown(a,i,n),n,ex,b[i], 1);
    }

  }
  else {

    int i;
    for (i=b.Size-1;i>=0;i--) {
      omuln(bitpown(a,i,n),n,ex,b[i],-1);
    }

    ~ ex[0];

  }

}

Integer shornuf(Integer Number) {

  Integer Factor;
  Integer f1,f2;
  Integer A;
  Integer Cv;
  Integer Xv;

  int size = ilog2(Number);

  cout << "Bits : " << size << "\n";

  QuantumRegister F(size);
  Integer Fv;

  QuantumRegister X(size);
  Integer Pv;

  int attempt = 1;
  int found = 0;
  while ((attempt < 8) && (! found)) {

    cout << "Attempt # " << attempt << "\n" ;

    do {
      
      A =  (rand(size) % Number) | 3 ;
    
    } while (gcd(Number,A) != 1) ;

    cout << "Random A : " << A.Value() << "\n";

    X.Set(0);

    X.Mix();
	    
    F.Set(0);

    expn(A,Number,X,F, 1); // F = a^X mod n

    Fv = F.Sample();

    Pv = X.Period();

    f1 = gcd((powmod(A,Pv/2,Number) + 1) % Number , Number) ;
    f2 = gcd((powmod(A,Pv/2,Number) + Number - 1) % Number , Number) ;

    Xv = X.Sample();

    Cv = powmod(A,Xv,Number);

    qdd_assert(Cv == Fv, "Mismatch");

    f1 = gcd(Number,f1);
    f2 = gcd(Number,f2);

    Factor = (f1 > f2 ? f1 : f2);

    found = (Factor > 1) && (Factor < Number);

    attempt++;

  }

  return(Factor);

}

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  Main

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

int main() {

  qdd_init(16,45); // Initialize for 16 bits and ~45 MEG of main memory

  try {

      Integer Number = 0;
      int value;

      while ((Number <= 8) || (Number > 0xFFFF)) {

	cout << "Enter Number to be factored (8 <= Number <= " << 0xFFFF << ") : " ;

	cin >> value;

        Number = value;

      }

      cout << "Factoring : " << Number.Value() << "\n" ;

      time_t t0;

      t0 = time(0);

      Integer Factor = shornuf(Number);

      Integer Alt = (Number / Factor);

      cout << "Done : " << Number.Value() << " = " << Factor.Value() << " * " << Alt.Value() << "\n";

      cout << "(time = " << (time(0) - t0) << " sec.)\n" ;

  }
  catch (QuantumException e) {
    qdd_error(e.msg);
  }

  return(0);
  
}


